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ABSTRACT 

In 2009, the discovery of two planets orbiting the evolved binary star system HW Virginis 
was announced, based on systematic variations in the timing of eclipses between the two stars. 
The planets invoked in that work were significantly more massive than Jupiter, and moved on 
orbits that were mutually crossing - an architecture which suggests that mutual encounters 
and strong gravitational interactions are almost guaranteed. 

In this work, we perform a highly detailed analysis of the proposed HW Vir planetary 
system. First, we consider the dynamical stability of the system as proposed in the discovery 
work. Through a mapping process involving 91,125 individual simulations, we find that the 
system is so unstable that the planets proposed simply cannot exist, due to mean lifetimes of 
less than a thousand years across the whole parameter space. 

We then present a detailed re-analysis of the observational data on HW Vir, deriving a 
new orbital solution that provides a very good fit to the observational data. Our new analysis 
yields a system with planets more widely spaced, and of lower mass, than that proposed in 
the discovery work, and yields a significantly greater (and more realistic) estimate of the 
uncertainty in the orbit of the outermost body. Despite this, a detailed dynamical analysis of 
this new solution similarly reveals that it also requires the planets to move on orbits that are 
simply not dynamically feasible. 

Our results imply that some mechanism other than the influence of planetary companions 
must be the principal cause of the observed eclipse timing variations for HW Vir If the sys- 
tem does host exoplanets, they must move on orbits differing greatly from those previously 
proposed. Our results illustrate the critical importance of performing dynamical analyses as a 
part of the discovery process for multiple-planet exoplanetary systems. 

Key words: binaries: close, binaries: ecUpsing, stars: individual: HW Vir, planetary systems, 
methods: N-body simulations 



1 INTRODUCTION 



Since the discovery of the first planets around other stars ( |Wol-| 
[szczan & Frail|1992| [Mayor & Queloz|1995| l, the search for exo- 
planets has blossomed to become one of the most exciting fields of 
modem astronomical research. The great majority of the hundreds 
of exoplanets that have been discovered over the past two decades 
have been found orbiting Sun-like stars by dedicated international 
radial velocity programs. Among these programmes are HARPS 
(the High Accuracy Radial velocity Planetary Search project, e.g. 
|Pepe et al.|2004i|Udry et al.|2007[payoretal.|2009j , AAPS (the 
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Anglo-Australian Planet Search, e.g. [Tinney et al.|2001[|2011[|Wit 



tenmyer et al. |2012bll, California (e.g. [Howard et al.|2010 Wright 



20n), Lick^amegie (e.g. 'Rivera et al.'"2010' 'Anglada- 



Escude et al.|2012}, and Texas (e.g. Endl et al. 2006; Robertson 



et al.[2012at . The other main method for exoplanet detection is the 



transit technique, which searches for the small dips in the bright- 
ness of stars that result from the transit of planets across them. 
Ground-based surveys such as WASP fHellier et al."201 1 /Smith et[ 
[al.[2012[ l and HAT ^Bakos et al.[2007l|Howard et al. 2012 ) have pi- 
oneered such observations, and resulted in the discovery of a num- 
ber of interesting planetary systems. In the coming years, such sur- 
veys using space-based observatories will revolutionise the search 
for exoplanets. Indeed, a rapidly growing contribution to the cata- 
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logue of known exoplanets comes from the Kepler spacecraft (e.g. 
IBorucki et al.|2011[, Welsh et al.|2012{ [Doyle et al. 2011 1, which 
will likely result in the number of known exoplanets growing by an 
order of magnitude in the coming years. 

In recent years, a number of new exoplanet discoveries have 
been announced featuring host stars that differ greatly from the 
Sun-like archetype that make up the bulk of detections. The most 
striking of these are the circumbinary planets, detected around 
eclipsing binary stars via the periodic variations in the timing of 
observed stellar eclipses. A number of these unusual systems fea- 
ture cataclysmic variable stars, interacting binary stars composed of 
a white dwarf primary and a Roche lobe filling M star secondary. 
Sharply-defined eclipses of the bright accretion spot, with periods 
of hours, can be timed with a precision of a few seconds. In these 
systems, the eclipse timings are fitted with a linear ephemeris, and 
the residuals (O — C) are found to display further, higher-order vari- 
ations. These variations can be attributed to the gravitational effects 
of distant orbiting bodies which tug on the eclipsing binary stars, 
causing the eclipses to appear slightly early or late. This light-travel 
time (LTT) effect can then be measured and used to infer the pres- 
ence of planetary-mass companions around these highly unusual 
stars. Some examples of circumbinary companions discovered in 



this manner include UZ For (Potter et al. 



2011b, NN Ser (Beuer- 



mann et al. 20 10)l, DP Leo jQian et al.|2010^, HU Aqr jSchwarz 



et al..2009,.Qian et al.|201I| i and SZ Her ( |Hinse et al.|20I2b[|L'ee 



et al.|2012| > 



The first circumbinary planets to be detected around hosts 
other than pulsars were those in the HW Vir system ( |Lee et aT] 
|2009[ l, which features a subdwarf primary, of spectral class B, and 
a red dwarf companion which display mutual eclipses with a pe- 
riod of around 2.8 hours jMenzies & Marang|1986| l. The detection 
of planets in this system was based on the timing of mutual eclipses 
between the central stars varying in a fashion that was best fit by 
including two sinusoidal timing variations. The first, attributed to a 
companion of mass M sin ! = 19.2 My„p, had aperiod of 15.8 years, 
and a semi-amplitude of 77 s, while the second, attributed to a com- 
panion of mass M sin ( = 8.5 Mjup, had a period of 9.1 years and 
semi-amplitude of 23 s. Whilst these semi-amplitudes might ap- 
pear relatively small, the precision with which the timing of mutual 
eclipses between the components of the HW Vir binary can be mea- 
sured means that such variations are relatively easy to detect. 

Over the last decade, a number of studies have shown that, for 
systems that are found to contain more than one planetary body, a 
detailed dynamical study is an important component of the planet 
discovery process that should not be overlooked (e.g. 'Stepinski et' 
ar]|2000[ [Gozdziewsk i & Maciejewski 2001 , Ferraz-Mello et al. 



2005 [ Laskar & Correia 20091. However, despite these pioneering 



works, the great majority of exoplanet discovery papers still fail to 
take account of the dynamical behaviour of the proposed systems. 
Fortunately, this situation is slowly changing, and recent discovery 
papers such as Robertson et al.|p012al, Wittenmyer et al. (2012b) 



and |Robertson et 3^(20 12b I have shown how studying the dynam- 
ical interaction between the proposed planets can provide signifi- 
cant additional constraints on the plausible orbits allowed for those 
planets. Such studies can even reveal systems in which the observed 
signal cannot be explained by the presence of planetary compan- 
ions. The planets proposed to orbit HU Aqr are one such case, with 



a number of studies (e.g. Homer et al.||2011 


Funk et al.|2011| 


Hinse et al.||20I2a| | Wittenmyer et al.||20I2a 


Gozdziewski et al.| 


20I2|) showing that the orbital architectures allowed by the observa- 



tions are dynamically unstable on astronomically short timescales. 
In other words, whilst it is clear that the observed signal is truly 



Parameter 


HW Vir b 


HW Virc 


Unit 


M sin (■ 


0.00809 ± 0.00040 


0.01836 ± 0.000031 




a sin i 


3.62 ± 0.52 


5.30 ± 0.23 


AU 


e 


0.31 ±0.15 


0.46 ± 0.05 




CO 


60.6 ±7.1 


90.8 ± 2.8 


deg 


T 


2,449,840 ± 63 


2,454,500 ± 39 


HJD 


P 


3316 ±80 


5786 ±51 


days 



Table 1. Parameters for the two planetary bodies proposed in the HW Vir 
system, taken from |Lee et al.|j2009^ (their Table 7). 



there, it seems highly unlikely that it is solely the result of orbiting 
planets. 

Given these recent studies, it is clearly important to consider 
whether known multiple-planet exoplanetary systems are truly 
what they seem to be. In this work, we present a re-analysis of 
the 2-planet system proposed around the eclipsing binary HW Vir. 
In section 2, we briefly review the HW Vir planetary system as 
proposed by |Lee et al.| ( |2009| . In section 3, a detailed dynamical 
analysis of the planets proposed in that work is perforined. In sec- 
tion 4, we present a re-analysis of the observations of HW Vir that 
led to the announcements of the exoplanets, obtaining a new orbital 
solution for those planets which is dynamically tested in section 5. 
Finally, in section 6, we present a discussion of our work, and draw 
conclusions based on the results herein. 

2 THE HW VIR PLANETARY SYSTEM 

The HW Vir system consists of a subdwarf B primary and an M6- 
7 main sequence secondary. The system eclipses with a period of 
2.8 hours ( [Menzies & Marang]|1986^ , and the stars have masses 
of Ml = 0.48 Mg and M2 = 0.14 Mq ( [Wood & Saffer|p^99l ). 
Changes in the orbital period of the eclipsing binary were first noted 
by Kilken ny et al.| ( [T994| ); further observations led other authors 
to suggest that the period changes were due to LTT effects aris- 
ing from an orbiting substellar comp anion ( [Kilkenny et al.|2003[ 



Ibanoglu et al. 



2004 



Lee et al. 



2009 1 obtained a further 8 years of 



photometric observations of HW Vir. Those data, in combination 
with the previously published eclipse timings spanning 24 years, 
indicated that the period changes consisted of a quadratic trend 
plus two sinusoidal variations with periods of 15.8 and 9.1 years. 
[Lee et al^ ( [2009^ examined alternative explanations for the cyclical 
changes, ruling out apsidal motion and magnetic period modula- 
tion via the Applegate mechanism ( [Applegate||1992^ . They con- 
cluded that the most plausible cause of the observed cyclic pe- 
riod changes is the light-travel time effect induced by two com- 
panions with masses 19.2 and 8.5 My,,^. The parameters of their 
fit can be found in Table[T| Formally, the planets are referred to as 
HW Vir (AB)b and HW Vir (AB)c, but for clarity, we refer to the 
planets as HW Vir b and HW Vir c. 

A first look at the fitted parameters for the two proposed 
planets reveals an alarming result: the planets are both massive, 
in the regime that borders gas giants and brown dwarfs, and oc- 
cupy highly eccentric, mutually crossing orbits with separations 
that guarantee close encounters - generally a surefire recipe for 
dynamical instability (as seen for the proposed planetary system 
around HU Aqr (e.g. iHorner et al.|201 l[[Hinse et al.|20T2a||Mt^ 
[tenmyer et al. [20 12a| [Gozdziewski et al.|2012^ . 

3 A DYNAMICAL SEARCH FOR STABLE ORBITS 



The work by [Lee et al.[^2009) derived relatively high masses (8.5 
and 19.2 Mj^p) and orbital eccentricities (0.31 and 0.46), and so 



© 2012 RAS, MNRAS 000,[T]-?? 



A Dynamical Analysis of the Proposed Circumbinary HW Virginis Planetary System 3 



significant mutual gravitational interactions are expected. To assess 
the dynamical stability of the proposed planets in the HW Vir sys- 
tem, we performed a large number of simulations of the planetary 
system, following a successful strategy used on a number of pre 
vious studies (e.g. [M arshall et al.|2010| [Homer et al.|2011|[Mt 



tenmyer et al.||2012a[ [Wittenmyer e t al. 2012b; Robertson et 



al. 



2012a|b[ Homer et al.|2012 1. We used the Hybrid integrator within 



the n-body dynamics package MERCURY ( [Chambers & Miglior-| 
[ini|1997| [Chambers |1999| l, and followed the evolution of the two 
giant planets proposed by Lee et al.| ( [2009[ > for a period of 100 mil- 
lion years. In order to examine the full range of allowed orbital 
solutions, we composed a grid of plausible architectures for the 
HW Vir planetary system, each of which tested a unique combi- 
nation of the system's orbital elements, spanning the ± 3-CT range 
in the observed orbital parameters. Following our earlier work, the 
initial orbit of HW Vir c (the planet with the best constrained or- 
bit in [Lee et al.[j200"9> was held fixed at its nominal best-fit values 
(i.e. a = 5.3 AU, e = 0.46 etc.). The initial orbit of HW Vir b was 
varied systematically such that the full 3-CT error ellipse in semi- 
major axis, eccentricity, longitude of periastron and mean anomaly 
were sampled. In our earlier work, we have found that the two main 
drivers of stability or instability were the orbital semi-major axis 
and eccentricity (e.g. [Homer et al.|20Il) , and so we sampled the 
3-a region of these parameters in the most detail. 

In total, 45 distinct values of initial semi-major axis were 
tested for the orbit of HW Vir b, equally distributed across the 
full ± 3-CT range of allowed values. For each of these unique semi- 
major axes, 45 distinct eccentricities were tested, evenly distributed 
across the possible range of allowed values (i.e. between eccentrici- 
ties of 0.00 and 0.76). For each of the 2025 a — e pairs tested in this 
way, fifteen unique values of CO were tested (again evenly spread 
across the ± 3-CT range, whilst for each of the a — e — CO values, three 
unique values of mean anomaly were considered. In total, therefore, 
we considered 91,125 unique orbital configurations for HW Vir b, 
spread in a 45 x 45 x 15 x 3 grid ina — e — CO — M space. In each 
of our simulations, the masses of the planets were set to their min- 
imal Msini values, in order to maximise the potential stability of 
their orbits. The orbital evolution of the planets was followed for 
a period of 100 million years, or until one of the planets was ei- 
ther ejected (defined by that planet reaching a barycentric distance 
of 20 AU), a collision between the planets occurred, or one of the 
planets collided with the central stars. If such a collision/ejection 
event occurred, the time at which it happened was recorded. 

In this way, the lifetime of each of the unique systems was 
determined. This, in turn, allowed us to construct a map of the dy- 
namical stability of the system, which can be seen in Fig[6] As can 
be seen in that figure, none of the orbital solutions tested were dy- 
namically stable, with few a-e locations displaying mean lifetimes 
longer than 1,000 years. 

Remarkably, we find that the proposed orbits for the HW Vir 
planetary system are even less dynamically stable than those pro- 
posed for the now discredited planetary system around HU Aqr 
Homer et al.|20Ilt|Wittenmyeret al.|2012a[|Hinse et al.|20I2a| 



[Homer et al.[20I2 Gozdziewski et al.|2012 1. Simply put, our result 



' The four yellow/orange/red hotspots in that plot between a ~ 3.3 and 
a ~ 4.2 AU ai'e the result of four unusually stable runs, with lifetimes of 
120 kyr (a = 3.31 AU.e = 0.12), 250 kyr {a = 3.74 AU ,e = 0.22), 56 kyr 
(a = 4.02 AC/, e = 0.19) and 49 kyr (a = 4.24 AJ7,e = 0.05). Such "long- 
live" outliers are not unexpected, given the chaotic nature of dynamical in- 
teractions, but given the typically very short lifetimes observed can signifi- 
cantly alter the mean lifetime in a given bin. 



proves conclusively that, if there are planets in the HW Vir system, 
they must move on orbits dramatically different to those proposed 
by [Lee et al.| ( [2009} . This instability is not particularly surprising, 
given the high orbital eccentricity of planet c, which essentially en- 
sures that the two planets are on orbits that intersect one another, 
irrespective of the initial orbit of planet b. Given that the two plan- 
ets are not trapped within mutual mean-motion resonance (MMR), 
such an orbital architecture essentially guarantees that they will ex- 
perience strong close encounters within a very short period of time, 
ensuring the system's instability. 

4 ECLIPSE TIMING DATA ANALYSIS AND LTT MODEL 

Given the extreme instability exhibited by the planets proposed by 
[Lee et alT| p009), it seems reasonable to ask whether a re-analysis 
of the observational data will yield significantly different (and more 
reasonable) orbits for the planets in question. We therefore chose to 
re-analyse the observational data, following a similar methodology 
as applied in an earlier study of HU Aqr ( [Hinse et al.|20I2a) . 

At the basis of our analysis we use the combined mid-eclipse 
timing data set compiled by [Lee et al.| ( [2009} , including the times 
of secondary eclipses. The timing data used in [Lee et aL| p009[ > 
were recorded in the UTC (Coordinated Universal Time) time stan- 
dard, which is known to be non-uniform ( Bastian|2000[[Guinan &[ 
,Ribas 200 1 ). To eliminate timing variations introduced by acceler- 
ated motion within the Solar System, we therefore transformecj^the 
HJD (Heliocentric Julian Date) timing records in UTC time stan- 
dard into Barycentric Julian Dates (BID) within the Barycentric 
Dynamical Time (TDB) standard (Eastman et al .|20T0^ . A total of 
258 timing measurements were used spanning 24 years from Jan- 
uary 1984 (HJD 2 445 730.6) to May 2008 (HJD 2 454 607.1). We 
assigned 1-ct timing uncertainties to each data point by following 
the same approach as outlined in |Lee et ah] ( [2009) 1. 

For an idealised, unperturbed and isolated binary system, 
the linear ephemeris of future/past mid-eclipse (usually primary) 
events can be computed from 



(1) 



where E denotes the (independent) ephemeris cycle number, Tq is 
the reference epoch, and Pq measures the eclipsing period (~ 2.8 
hrs) of HW Vir. A linear regression performed on the 258 recorded 
light-curves allows Pq to be determined with high precision. In this 
work, we chose to place the reference epoch close to the middle 
of the observing baseline to avoid parameter correlation between 
T() and during the fitting process. In the following we briefly 
outline the LTT model as used in this work. 

4.1 Analytic LTT model 

The model adopted in this work is similar to that described in [Hinse[ 
[et al.| ( (2012a^ , and is based on the original formulation of a sin- 
gle light-travel time orbit introduced by ,Irwin et'ar] ( [T952| >. In this 
model the two components of the binary system are assumed to 
represent one single object with a total mass equal to the sum of 
the masses of the two stars. This point mass is then placed at the 
original binary barycentre. If a circumbinary companion exist, then 
the combined binary mass follows an orbit around the total system 
barycentre. The eclipses are then given by Eq.[T] This defines the 
LTT orbit of the binary. The underlying reference system has its 
origin at the total centre of mass. 



http://astroutils.astronomy.ohio-state.edu/time 
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Following [Irwin et al.| ^1952[ l, if the observed mid-eclipse 
times exhibit a sinusoidal-like variation (due to one or more un- 
seen companion(s)), then the quantity — C defines the light-travel 
time effect and is given by 



ephemeris of calculated times of mid-eclipses then takes the form 



{O^C){E) = To{E)-Tc{E)-. 



2 



(2) 



where Tq denotes the measured time of an observed mid-eclipse, 
and Tc is the computed time of that mid-eclipse based on a linear 
ephemeris. We note that Xi +X2 is the combined LTT effect from 
two separate two-body LTT orbits. The quantity X; is given by the 
following expression for each companion ( |Irwin et al.|1952| l: 



l-e 



h.i 



+ eh J COS fhj 



sin(/i./ + co^,,,) + ehj sinah i 



(3) 



where K^j = ai, isinl^j/c is the semi-amplitude of the light-time 
effect (in the — C diagram) with c measuring the speed of light 
and Ih i is the line-of-sight inclination of the LTT orbit relative 
to the sky plane, ei,j the orbital eccentricity, , the true longi- 
tude and CO/, ,- the argument of pericenter of the LTT orbit. The 
5 model parameters for a single LTT orbit are given by the set 
(ahjSinIi, i,ei,j,Oii, j,Ti,j,Ph,i). The time of pericentre passage T/, ,- 
and orbital period P;, , are introduced through the expression of 
the true longitude as a time-like variable via the mean anomaly 
M = ni,i(To — Ti,i), with nj,i = 2tz/Pi,j denoting the mean mo- 
tion of the combined binary in its LTT orbit. Computing the true 
anomaly as a function of time (or cycle number) requires the solu- 
tion of Kepler's equation. We direct the interested reader to |Hinse| 
|et al.| l |2012"al > for further details. 

In Eq. [3] the origin of the coordinate system is placed at the 
centre of the LTT orbit (see e.g. |Irwin et al.|1952y . A more natural 
choice (from a dynamical point of view) would be to use the system 
centre of mass as the origin of the coordinate system. However, 
the derived Keplerian elements are identical in the two coordinate 
systems (e.g. [Hinse et al.|201 2b). Finally, we note that our model 
does not include mutual gravitational interactions. We also only 
consider the combination of two LTT orbits from two circumbinary 
companions. 

From first principles, some similarities exist between the LTT 
orbit and the orbit of the circumbinary companion. First, the ec- 
centricities (e;, = e,) and orbital periods (P;, = P,) are the same. 
Second, the arguments of pericenter are 180° apart from one an- 
other (CO; = 180° — (Ohj). Third, the times of pericenter passage are 
also identical (T/, j = 7]). 

Information of the mass of the unseen companion can be ob- 
tained from the mass function given by 

4Tl^{abjsmlbj)^ 4Tt^{Khjcf (m/sin/,)^ 



GPli 



(mh + mi) 



1,2 



(4) 

The least-squares fitting process provides a measure for A"/, ,■ and 
Ph i, and hence the minimum mass of the companion can be found 
from numerical iteration. In the non-inertial astrocentric reference 
frame, with the combined binary mass at rest, the companion's 
semi-major axis relative to the binary is then calculated using Ke- 
pler's third law. 

In [Lee et al.[j2009| l the authors also accounted for additional 
period variations due to mass transfer and/or magnetic interactions 
between the two binary components. These variations usually occur 
on longer time scales compared to orbital period variations due to 
unseen companions. Following [Hilditch] ( [2001^ , the corresponding 



Tc = TQ+PQE + ^E^+Xi, 1=1,2 



(5) 



where [3 is an additional free model parameter and accounts for a 
secular modulation of the mid-eclipse times resulting from interac- 
tions between the binary components. Assuming the timing data of 
HW Vir are best described by a two-companion system, and to be 
consistent with [Lee et al.[ l [2009t , we have used Eq.[5]as our model 
which consists of 13 parameters. 

5 METHODOLOGY AND RESULTS FROM 
-PARAMETER SEARCH 

To find a stable orbital configuration of the two proposed circumbi- 
nary companions, we carried out an extensive search for a best-fit 
in parameter space. The analysis, methodology and technique 
follow the same approach as outlined in [Hinse et al.[ ( [2012a^ . Here 
we briefly repeat the most important elements in our analysis. 

We used the Levenberg-Marquardt least-square minimisation 
algorithm as implemented i n the IDlj^based software package 
MPFIlj^ ( Markwardt 2009 1. The goodness-of-fit statistic of 
each fit was evaluated from the weighted sum of squared errors. 
In this work, we use the reduced chi-square statistic j^. which takes 
into account the number of data points and the number of freely 
varying model parameters. 

We seeded 28,201 initial guesses within a Monte Carlo experi- 
ment. Each guess was allowed a maximum of 500 iterations before 
termination, with all 13 model parameters (including the secular 
term) kept freely varying. Converged solutions with ^ 10 were 
accepted with the initial guess and final fitting parameters recorded 
to a file. After each converged iteration we also solved the mass 
function (Eq.|4]l for the companion's minimum mass and calculated 
the semi-major axis relative to the system barycentre from Kepler's 
third law. 

Initial guesses of the model parameters were chosen at random 
following either a uniform or normal distribution. For example, ini- 
tial orbital eccentricities were drawn from a uniform distribution 
within the interval e e [0.0,0.8]. Our initial guesses for the orbital 
periods were guided by a Lomb-Scargle (LS) discrete Fourier trans- 
formation analysis on the complete timing data set. For the LS anal- 
ysis we used the PERIOD04 software package ( [Lenz & Breger 
[2005[ l capable of analysing unevenly sampled data sets. Fig. [T 
shows the normalised LS power-spectrum. The LS algorithm found 
two significant periods with frequencies /i ~ 1.4 x 10^"* cycles/day 
and /2 ~ 2.1 X 10^"* cycles/day. These frequencies correspond to 
periods of 7397 and 4672 days, respectively. Hence the short-period 
variation is covered more than twice during the observing interval. 
Due to a lower amplitude, it contains less power within the data set. 
Our random initial guesses for the companion periods were then 
drawn from a Gaussian distribution centred at these periods with 
standard deviation of ± 5 years. We call this approach a "quasi- 
global" search of the underlying x^-parameter space. 

5.1 Results - finding best fit and confidence levels 

Our best fit model resulted in a Xr = 0.943 and is shown in Fig. [2] 
along with the LTT signal due to the inner and outer companion 
and the secular term. The corresponding root-mean-square (RMS) 

^ The acronym IDL stands for Interactive Data Language and is a 
trademark of ITT Visual Information Solutions. For further details see 
Jittp : //www. ittvis . com/Product Service 5 /IDL . aspx^ 
http://purl.com/net/mptit 
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scatter of data around the best fit is 8.7 seconds, which is close 
to the RMS scatter reported in |Lee et af] ( |2009 l. The fitted model 
elements and derived quantities of our best fit are shown in Ta- 
ble |2] Compared to the system in |Lee et af] ( |2009^ we note that 
we now obtain a lower eccentricity for the inner companion and a 
larger eccentricity for the outer companion. Furthermore, our two- 
companion system has also slightly expanded, with larger semi- 
major axes (and therefore longer orbital periods) for both compan- 
ions compared to |Lee et aL| l |2009| l. 

The next question to ask is how reliable or significant our 
best-fit solution is in a statistical sense. Assuming that the errors 
are normally distributed one can establish confidence levels for a 
multi-parameter fit jBevington & Robinson|[r992^ . We therefore 
carried out detailed two-dimensional parameter scans covering a 
large range around the best-fit value in order to study the x^-space 
topology in more detail. In particular, we explored relevant model 
parameter combinations including T{),P{) and (3. 

In all our experiments we allowed the remaining model param- 
eters to vary freely while fixing the two parameters of interest in the 
considered parameter range ( B evington & Robinson|[T992| |Press| 
|et al.|I992) . Assuming parameter errors are normally distributed, 
our I-, 2- and 3-CT level curves provide the 68.3, 95.4 and 99.7% 
confidence levels relative to our best-fit, respectively. In Fig.|3]we 
show a selection of our two-dimensional parameter scan consider- 
ing various model parameters. 

Ideally, one would aim to work with parameters with little 
correlation between the two parameters. The lower-right panel in 
Fig.|3]shows the relationship between To and P^. The near circular 
shape of the level curves reveals that little correlation between the 
two parameters exists. This is most likely explained by our choice 
to locate the reference epoch in the middle of the dataset. The re- 
maining panels in Fig. [3] show some correlations between the pa- 
rameters. However, we have some indication of an unconstrained 
outer orbital period in the lower-left panel of Fig.|3] We show the 
location of orbital mean motion resonances in the {P\,P2) plane. 
Our best fit is located close to the 2: 1 mean motion resonance. 

To demonstrate that the outer period is unconstrained we have 
generated /^-parameter scans in the {P2,e2)-plane as shown in 
Fig.|4] Our best fit model is shown in the left-most panel of Fig.|4] 
along with the l-,2- and 3-0 confidence levels. It is readily evident 
that the l-O confidence level does not simply surround our best-fit 
model in a confined or ellipsoidal manner. We rather observe that 
all three level curves are significantly stretched towards solutions 
featuring longer orbital periods for the outer companion. This is 
demonstrated in the middle and right panels of Fig. |4] Any longer 
orbital period for the outer companion therefore results in a with 
similar statistical significance as our best-fit model. 

In addition, we studied the (asinz,f )-parameter plane for the 
outer companion, as shown in the lower-middle panel of Fig. |3] 
Here, we also observe that asini is unconstrained. The uncon- 
strained nature of the outer companion's asini and orbital period 
has a dramatic effect on the derived minimum mass and the corre- 
sponding error bounds. 

5.2 Results - parameter errors 

When applying the LM algorithm formal parameter errors are ob- 
tained from the best-fit covariance matrix. However, in our study, 
we sometimes encounter situations where some of the matrix el- 
ements become zero, or have singular values. However, at others 
times, the error matrix is returned with non-zero elements. In those 
cases we have observed that the outer orbital period is often better 



determined than that of the inner companion. In the case of HW Vir, 
such solutions are clearly incongruous, given the relatively poor 
orbital characterisation of the outer body compared to that of the 
innermost. 

Being suspicious about the formal covariance errors, we have 
resorted to two other methods to determine parameter errors. First, 
we attempted to determine errors by the use of the bootstrap method 
( [Press et al.||1992"l l. However, we found that the resulting error 
ranges are comparable to the formal errors extracted from the best- 
fit covariance matrix. Furthermore, the bootstrap error ranges were 
clearly incompatible with the l-o error "ellipses" discussed above. 
For example, from the top-left panel in Fig. [3] we estimate the l-c 
error on the inner orbital period to be on the order of 50-100 days. 
In contrast, the errors for the inner orbital period obtained from our 
bootstrap method were of order just a few days. For this reason, we 
consider the bootstrap method to have failed, and it has therefore 
not been investigated further in this study. However, it is interest- 
ing to speculate on the possibility of dealing with a dataset which 
is characterised by "clumps of data", as seen in Fig.|2] When gen- 
erating random (with replacement) bootstrap ensembles, there is 
the possibility that only a small variation is being introduced for 
each random draw, as a result of the dumpiness of the underlying 
dataset. That dumpiness might be mitigated for, and the bootstrap 
method rendered still viable for the establishment of reliable errors, 
by enlarging the number of bootstrap ensembles to compensate for 
the lack of variation within single bootstrap data sets. One other 
possibility would be to replace clumps of data by a single data point 
reflecting the average of the clump. However, we have instead in- 
voked a different approach. 

Having located a best-fit minimum, we again seeded a large 
number of initial guesses around the best-fit parameters. This time, 
we considered only a relatively narrow range around the best-fit 
values (e.g. as shown in Fig. [3] and Fig. |4] This ensured that the 
LM algorithm would iterate towards our best-fit model depend- 
ing on the underlying topology and inter-parameter correlations. 
The initial parameter guesses were randomly drawn from a uniform 
distribution within a given parameter interval. We then iterated to- 
wards a best-fit value using LM, and recorded the best-fit parame- 
ters along with the corresponding Xr- 

Generating a large number of guesses enabled us to establish 
statistics on the final best-fit parameters with x^ within I-, 2-, and 
3-CT confidence levels. We therefore performed a Monte Carlo ex- 
periment that considered several tens of thousands of guesses. To 
establish l-a error bounds (assuming a normal distribution for each 
parameter), we then considered only those models that yielded Xr 
within the l-a confidence limits (inner level curves), as shown in 
Fig. [3] and Fig.|4] The error for a given parameter is then obtained 
from the mean and standard deviation, and listed in Table[2] In or- 
der to test our assumption of normally distributed errors, we plotted 
histograms for the various model parameters in Fig.|5] For each his- 
togram distribution, we fitted a Gaussian and established the corre- 
sponding mean and standard deviation. While some parameters fol- 
low a Gaussian distribution (for example the outer companion's ec- 
centricity), other parameters show no clear sign of "Gaussian tails". 
This is especially true for the orbital period of the outer compan- 
ion. Following two independent paths of analysis, we have demon- 
strated that the outer period is unconstrained, based on the present 
dataset, and should be regarded with some caution. However, we 
also point out a short-coming of our method of determining random 
parameter errors. The parameter estimates depend on the proximity 
of the starting parameters to the best-fit parameter. In principle, our 
method of error determination assumes that the best-fit model pa- 
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rameters are well-determined in terms of well-established closed- 
loop confidence levels around the best-fit parameters. The true ran- 
dom parameter error distribution for the two ill-constrained param- 
eters might turn out differently. To put stronger constraints on the 
model parameters is clearly only possible by augmenting the ex- 
isting timing data through a program of continuous monitoring of 
HW Vir over the coming years (Pribulla et al.|2012||Konacki et al.^ 
|2012| l. As more data is gathered, the confidence levels in Fig.|4]will 
eventually narrow down. 



6 THE P COEFFICIENT AND ANGULAR MOMENTUM 
LOSS 

One of the features of our best-fit orbital solution is that it results in 
a relatively large p coefficient, which can be related to a change in 
the period of the binary resulting from additional, non-planetary ef- 
fects. Potential causes of such a period change include mass trans- 
fer, loss of angular momentum, magnetic interactions between the 
two binary components and/or perturbations from a third body on 
a distant and unconstrained orbit. In this study, the p factor repre- 
sents a constant binary period change (see |Hilditch|200l"| p. 171) 
with a linear rate of dP/dt = —9.57 x 10^^ days/yr, which is about 
15 percent larger than the value reported in |Lee et al.H2009l l. We 
retained the p coefficient in our model in order that our treatment 
be consistent with that detailed in |Lee et al.| ( |2009l l, such that our 
results might be directly compared to their work. 

|Lee et al!] t2009 ) examined a number of combinations of mod- 
els that incorporated a variety of potential causes for the observed 
period modulation. They found that the timing data is best de- 
scribed by two LTT and a quadratic term in the linear ephemeris 
model. In their work, |Lee et ah] |2009| carefully examined the 
contribution of period modulation by various astrophysical effects. 
They were able to provide arguments that rule out the operation of 
the Applegate mechanism, due to the lack of small-scale variations 
in the observed luminosities that would have an influence on the 
J2 oblateness coefficient of the magnetically active component. A 
change in J2 would, in turn, affect the binary period. 

Furthermore, [Lee et aT] ^2009^ reject the idea that the ob- 
served O — C variation could be the result of apsidal motion, based 
on the circular orbit of the HW Vir binary system. In addition, 
they estimated the secular period change of the HW Vir binary or- 
bit due to angular momentum loss through gravitational radiation 
and magnetic breaking. They found that the most likely explana- 
tion for the observed linear decrease in the binary period is that it 
is the result of angular momentum loss by magnetic stellar wind 
breaking in the secondary M-type component. From first princi- 
ples (e.g. [Brinkworth et al.~||2006[ l, the period change observed 
in this work corresponds to a angular momentum change of order 
dj /dt = —2.65 X 10^^ erg. This is approximately 15 percent larger 
than the value reported by |Lee et al.| (2009 1, but still well within the 
range where magnetic breaking is a reasonable astrophysical cause 
for period modulation. 

Finally, we note that, whilst this work was under review, 
|Beuermann et al.|j2012^ published a similar study based on new 
timing data, in which they also considered the influence of period 
changes due to additional companions. In their work, they obtained 
a markedly different orbital solution than those discussed in this 
work, one which they found to be dynamically stable on relatively 
long timescales. In light of their findings, it is interesting to note 
that they do not include a quadratic term in their linear ephemeris 
model. This could point at the possibility that the inclusion of a 



quadratic term is somehow linked to the instability of the best-fit 
system found in this work. 

[Beuermann et al.H2012[ l found stable orbits for a solution in- 
volving two circumbinary companions. However, despite this, we 
note that the two models share some qualitative characteristics. A 
careful examination of Fig. 2 in Beuerma nn et al.| ( |2()12) l reveals 
that the orbital period of the outer companion is unconstrained 
from a period analysis since the x^-contour curves are open towards 
longer outer orbital periods - a result mirrored in our current work. 
As we noted earlier, we were unable to place strict confidence lev- 
els on the best-fit outer companion's orbital period and semi-major 
axis. Although [Beuermann et al.]p012| > do find a range of stable 
scenarios featuring their outer companion, we note that they fixed 
the eccentricity of that companion's orbit to be near-circular, with 
period of 55 years. Such an assumption (i.e. fixing some orbital pa- 
rameters) is somewhat dangerous, since it can lead to the produc- 
tion of dynamically stable solutions that are not necessarily sup- 
ported by the observational data (e.g. [Homer et al.|2012^ . A more 
rigorous strategy would be to generate an ensemble of models with 
each model (all parameters freely varying) tested for orbital stabil- 
ity (using some criterion like non-crossing orbits or non-overlap of 
mean-motion resonances, etc.) resulting in a distribution of stable 
best-fit models. Using the new data set, a study of the distribution 
of the best-fit outer planet's eccentricity would be interesting. It is 
certainly possible that the new data set constrains this parameter 
sufficiently in order to validate their assumptions. Although the re- 
sults presented in |Beuermann et aL] ( |2012^ are clearly promising, 
it is definitely the case that more observations are needed before 
the true origin of the observed variation for HW Vir is established 
beyond doubt. 

The question of whether a period damping factor is truly nec- 
essary for the HW Vir system would require a statistically self- 
consistent re-examination of the complete data set taking account 
of a range of model scenarios. We refer the interested reader to 
[Gozdziewski et aL| ( [2012[ l, who recently carried out a detailed in- 
vestigation of the influence of the quadratic term for various sce- 
narios in their attempt to explain the timing data of the HU Aqr 
system. 

7 DYNAMICAL ANALYSIS OF THE BEST-FIT LTT 
MODEL 

Since a detailed re-analysis of the observational data on the HW Vir 
system yields a new orbital solution for the system, it is interesting 
to consider whether that new solution offers better prospects for dy- 
namical stability than that proposed in |Lee et al l ( 2009 1 . We there- 
fore repeated our earlier dynamical analysis using the new orbital 
solution. Once again, we held the initial orbit of planet HW Vir c 
fixed at the nominal best-fit solution, and ran an equivalent grid of 
unique dynamical simulations of the planetary system, varying the 
initial orbit of HW Vir b such that a total of 45 distinct values of a 
and e, 15 distinct values of CO and 3 values of M were tested, each 
distributed evenly as before across the ± 3-CJ range of allowed val- 
ues. As before, the two simulated planets were assigned the nomi- 
nal M sin i masses obtained from the orbital model. The results of 
our simulations can be seen in Fig.|7] 

As was the case for the original orbital solution proposed in 
|Lee et af] ( |2009[ l, and despite the significantly reduced uncertain- 
ties in the orbital elements for the resulting planets, very few of the 
tested planetary systems survived for more than 1 ,000 years (with 
just twenty six systems, 0.029% of the sample, surviving for more 
than 3,000 years, and just three systems surviving for more than 
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10,000 years). As was the case for the planetary system proposed 
to orbit the cataclysmic variable system HU Aqr (e.g. [Horner et al.| 
|2011[[Wittenmyer et al.|2012a^ , it seems almost certain that the pro- 
posed planets in the HW Vir system simply do not exist - at least on 
orbits resembling those that can be derived from the observational 
data. 



8 CONCLUSION AND DISCUSSION 

The presence of two planets orbiting the evolved binary star sys- 
tem HW Vir was proposed by |Lee et aL| j2009^ , on the basis of 
periodic variations in the timing of eclipses between the two stars. 
The planets proposed in that work were required to move on rela- 
tively eccentric orbits in order to explain the observed eclipse tim- 
ing variations, to such a degree that the orbit of the outer planet 
must cross that of the innermost. It is obvious that, when one ob- 
ject moves on an orbit that crosses that of another, the two will 
eventually encounter one another, unless they are protected from 
such close encounters by the influence of a mutual mean-motion 
resonance (e.g. [Homer et al.|2004a|b| l. Even objects protected by 
the influence of such resonances can be dynamically unstable, al- 
beit on typically longer timescales (e.g. [H omer & Lykawka|2010| 
[Homer, Miiller & Lykawka |2012[ [Horne r et al. 2012}. Since the 
two planets proposed by Lee et al!]l(2009 ) move on calculated or- 



bits that allow the them to experience close encounters and yet are 
definitely not protected from such encounters by the influence of 
mutual mean-motion resonance, it is clear that they are likely to be 
highly dynamically unstable. To test this hypothesis, we performed 
a suite of highly detailed dynamical simulations of the proposed 
planetary system to examine its dynamical stability as a function of 
the orbits of the proposed planets. We found the proposed system 
to be dynamically unstable on extremely short timescales, as was 
expected based on the proposed architecture for the system. 

Following our earlier work jWittenmyer et al.|2012a[ |Hinse[ 
[et al.[|2012a^ , we performed a highly detailed re-analysis of the 
observed data, in order to check whether such improved analy- 
sis would yield better constrained orbits that might offer better 
prospects for dynamical stability. Our analysis resulted in calcu- 
lated orbits for the candidate planets in the HW Vir system that have 
relatively small uncertainties. Once these orbits had been obtained, 
we performed a second suite of detailed dynamical simulations to 
ascertain the dynamical stability of the newly determined orbits. 
Following the same procedure as for the original orbits, we con- 
sidered the stability of all plausible architectures for the HW Vir 
system. Despite the increased precision of the newly determined 
orbits, we find that the planetary system proposed is dynamically 
unstable on timescales as short as a human lifetime. For that rea- 
son, we must conclude that the eclipse-timing variations observed 
in the HW Vir system are not solely down to the gravitational in- 
fluence of perturbing planets. Furthermore, if any planets do exist 
in that system, they must move on orbits dramatically different to 
those considered in this work. 

Our results highlight the importance of performing comple- 
mentary dynamical studies of any suspected multiple-exoplanet 
system - particularly in those cases where the derived planetary or- 
bits approach one another closely, are mutually crossing and/or de- 
rived companion masses are large. Following a similar strategy as 
applied to the proposed planetary system orbiting HU Aqr i Homer[ 



[et al.|2011[[Wittenmyer et al.|2012"a)[Gozdziewski et al.|2012^ , we 
have found that the proposed 2-planet system around HW Vir does 
not stand up to a detailed dynamical scmtiny. In this work we have 
shown that the outer companion's period (among other parameters) 



is heavily unconstrained by establishing confidence limits around 
our best-fit model. However, we also point out the fact that the 
two circumbinary companions have brown-dwarf masses. Hence, 
a more detailed n-body LTT model which takes account of mutual 
gravitational interactions might provide a better description of the 
problem. 

To further characterise the HW Vir system and constrain or- 
bital parameters we recommend further observations within a mon- 
itoring program as described in'Pribulla et al.[(2012[ l. In a recent 
work on HU Aqr, Gozdziewski et al. (2012]l pointed out the pos- 
sibility that different data sets obtained from different telescopes 
could introduce systematic errors resulting in a false-positive de- 
tection of a two-planet circumbinary system. 

Finally, we note that, whilst this paper was under referee, 
[Beuermann et al.| ( [2012^ independently published their own new 
study of the HW Vir system. Based on new observational timing 
data, those authors determined a new LTT model that appears to 
place the 2-planet system around HW Vir on orbits that display 
long-term dynamical stability. Based on the results presented in 
this work, we somewhat doubt their findings of a stable 2-planet 
system and question whether such a system is really supported by 
the new data set given that no strict confidence levels were found 
for the best-fit outer period. Since performing a full re-analysis of 
their newly compiled data, including dynamical mapping of their 
new architecture, would be a particularly time intensive process, 
we have chosen to postpone this task for a future study. 
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Figure 1. Power spectrum of HW Vir timing data set using BJD(TDB) times witli tlie O — C residuals measured in seconds. Two periods were found with 
= 7397 and 1 /f2 = 4672 days, corresponding to 20.3 and 12.8 years, respectively. S/N denotes the frequency signal-to-noise ratio, which are significantly 
larger than the spectrum's noise level. Normalisation was done by division of the maximum amplitude in each spectrum. The /2-frequency was determined 
from the residuals after subtracting the period I//1 from the original timing data set. Additional peaks in both panels represent 1-year alias frequencies due to 
the repeating annual observing cycle of HW Vir. 
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Figure 2. Best-fit two-Kepler LTT model with xf- = 0.943. The best-fit model parameters (including reference epoch) are shown in Table|2]RMS denotes the 
root-mean-square scatter around the best fit. The lower part of the figure shows the residuals between the best fit model and the observed timing data set. 
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Figure 3. Colour-coded Xr parameter scans of orbital pai'ameters with remaining parameters to vary freely. The best-fit parameter is indicated by a star-like 
symbol. Contour curves show the l,2,3o confidence level curves around our best-fit model. See electronic version for colour figures. 
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Figure 4. Same as Fig.jsjbut considering the (P2 1*2) -plane- Tlie middle and riglit panel shows the Xr topology for longer orbital periods of the outer companion. 
Based on the used data set, no firm confidence levels can be established around out best-fit value. 
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Figure 5. Histogram distribution of six model parameters as obtained from a Monte Carlo experiment. Only models with < 1 — o where considered to 
assess the 68% confidence levels for each parameter Solid curves show fitted normal distribution with mean and standard deviation indicated in each panel. 
However, we used the mean and standard deviation derived from the underlying dataset to derive our 1-a errors as quoted in Table|2] 



© 2012 RAS, MNRAS 000,[T]-?? 



12 J. Horner, T. C. Hinse, R.A. Wittenmyer, J. P. Marshall & C.G. Tinney 




2 2.5 3 3.5 4 4.5 5 

Semi -Major Axis (AU) 

Figure 6. The stability of the HW Vir planetary system as proposed by |Lee et al.H2009) , as a function of the semi-major axis, a, and eccentricity, e, of planet 
HW Vir b. The initial orbit of HW Vir c was the same in each integration, set to the nominal best fit orbit from that work. The mean lifetime of the planetary 
system (in log\Q(lif etime jyr)) at a given a — e co-ordinate is denoted by the colour of the plot. The lifetime at each a — e location is the mean value of 45 
separate integrations canied of orbits at that a — e position (testing a combination of 15 unique CO values, and 3 unique M values). The nominal best-fit orbit 
for HW Vir b is located within the open square, from which lines radiate showing the extend of the ± l-o errors on a and e. As can be seen, the orbits of the 
system are incredibly unstable, no matter what initial orbit is considered for HW Vir b. 
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Figure 7. The stability of the HW Vir planetary system, given the orbital solution derived in this work, as a function of the semi-major axis, a, and eccentricity, 
e, of planet HW Vir b. The initial orbit of HW Vir c was the same in each integration, set to the nominal best fit orbit as detailed in Table 2. The mean lifetime 
of the planetary system (shown as log\o(lifetime/yr) at a given a — e co-ordinate is denoted by the colour of the plot. The lifetime at each a — e location is 
the mean value of 45 separate integrations carried of orbits at that a — e position (testing a combination of 15 unique CO values, and 3 unique M values). The 
nominal best-fit orbit for HW Vir b is located within the open square, from which lines radiate showing the extend of the ± 1 -o errors on a and e. Once again, 
the orbits of the system are found to be incredibly unstable, no matter what initial orbit is considered for HW Vir b. The two red hotspots in that plot are the 
result of two unusually stable runs, with lifetimes of 33 kyr (a = 4.185 AU, e = 0.137) and 38 kyr (a = 4.15 AU, f = 0.1 1). Even these most extreme outliers 
are dynamically unstable on astronomically short timescales. 
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Parameter 




two-LTT 


Unit 






T2 {i = 2) 








0.943 




RMS 




8.665 


seconds 


P 


-1.529 - : 


10"'2±1.25-10"'3 


day/cycle" 


To 
Po 

eb,i (or(;i,2) 
to*,i 

n,i 

Phi (or Pi 2) 


2,450, 280.28596 ±2.3 - 10-5 
0.116719519±4.6- 10"' 
0.081 ±0.002 0.196±0.012 
0.17±0.02 0.61 ±0.02 
0.05 ±0.01 2.09 ±0.08 
2, 448, 880 ±57 2, 448, 629 ±42 
4021 ±64 7992±551(!) 


BJD(TDB) 
days 
AU 

rad. 
BJD(TT) 

days 


Kbj 


4.6- 10-''±1.3- 10 


1.13- Kr3±7.04- 10-5 


days 


nii sin/, 

Oj sin/, 

ei 

to, 

Ti 

Pi 


12±3 
4.26 ±0.05 
0.17±0.02 
(71 -0.05) ±0.01 
2,448, 880±57 
4021 ±64 


11±8 
6.8±0.3 
0.61 ±0.02 
(tc- 2.09) ±0.08 
2, 448, 629 ±42 
7992±551 (!) 


Mjup 
AU 

rad. 
BJD(TT) 
days 



Table 2. Best-fit parameters for the LTT orbits of HW Vir corresponding to Fig.|2] Subscripts 1,2 refer to the circumbinary companions with / = 1, the inner, 
and (' = 2, the outer, companions. RMS measures the root-mean-square scatter of the data around the best fit. 1-a uncertainties have been obtained as described 
in the text. The last five entries are quantities of the two companions in the astrocentric coordinate system. Note that our values for a and P are somewhat 
larger than those of |Lee et al.|j2009| : (if, = 3.62 and Uc = 5.30 AU, Pb = 3316 and P^ = 5786 days. 
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